#define MAX_COMPLEX_POINT_COUNT   6000000
#define MAX_N_SEGMENT			  (MAX_COMPLEX_POINT_COUNT/2)

class FarstConv
{
public:
	FarstConv()
	{
		//整个的数据输出长度
		int outLen	= sizeof(complex) * MAX_COMPLEX_POINT_COUNT;
		m_fftIn		= (complex *)malloc(outLen);

		//冲击函数fft缓存
		int hLen	= sizeof(complex) * MAX_N_SEGMENT;
		m_fftH		= (complex *)malloc(hLen);
		m_fftS		= (complex *)malloc(hLen);
		m_fftM		= (complex *)malloc(hLen);
		m_endC		= (complex *)malloc(hLen);
		m_hC		= (complex *)malloc(hLen);
	}

	/*
	N =	 原始数据序列长度
	M =  冲击函数序列长度
	L =  期望的滑窗分段长度
	*/
	void		doConv(complex* x, double* h, int N, int M, int L, complex* out, int & outLen)
	{
		// 每段总长度（FFT 长度)
		int N_segment	= L + M - 1;					
		N_segment		= nextpow2(N_segment);

		//修正滑窗长度
		L				= N_segment - M + 1;

		//在前面M-1个点进行补充0
		memset(&m_fftIn[0],			0,		sizeof(complex) * (M - 1));
		memcpy(&m_fftIn[M - 1],		&x[0],	sizeof(complex) * N);

		//补充0后数据点个数
		int AL			= (M - 1) + N;

		//分段个数
		int num_blocks	= ceil((double)(AL) / double(L));

		//冲击函数求fft（用复数fft以免出错）
		memset(&m_hC[0], 0, N_segment*sizeof(complex));
		

		{
			for (size_t i = 0; i < M; i++)
			{
				m_hC[i][0] = h[i];
			}
		}


		doCfft(&m_hC[0], N_segment, &m_fftH[0]);
		{
			outLen = 0;
			for (size_t i = 0; i < num_blocks; i++)
			{
				//以L进行滑窗，每段取N_segment个数据
				int startI = i* L;
				int endI = startI + N_segment - 1;

				if (endI >= AL) //如果结果的索引大于长度
				{
					//数据填充
					int lastLen = AL - startI;
					memcpy(&m_endC[0], m_fftIn[startI], sizeof(complex) * lastLen);

					//补0
					memset(&m_endC[lastLen], 0, sizeof(complex) * (N_segment - lastLen));
					doCfft(&m_endC[0], N_segment, &m_fftS[0]);
				}
				else
				{
					doCfft(&m_fftIn[startI], N_segment, &m_fftS[0]);
				}

				//做频域卷积
				mult_complex(&m_fftH[0], &m_fftS[0], N_segment, &m_fftM[0]);
				doIfft(&m_fftM[0], N_segment, m_endC);

				// 第三大步， 将每段的卷积结果，丢弃前M - 1个点，保留L个点，进行结果数据输出
				// 步骤5：提取有效数据（丢弃前M - 1个点）
				//int NLost				= N_segment - (M - 1);
				int bgIndex = L*i;
				outLen = outLen + L;
				memcpy(&out[bgIndex], m_endC[M - 1], sizeof(complex)*L);
			}
		}
	}

	int nextpow2(int N)
	{
		if (N <= 1) return 1;
		double log2x = log2(N);
		return pow(2, ceil(log2x));
	}

public:

	//做复数的fft
	void		doCfft(complex* x, int len, complex* out)
	{
		memset(&out[0], 0, sizeof(complex)* len);
		dft_1d(len, x, out, FFTW_FORWARD/*正向fft*/, FFTW_ESTIMATE);
	}

	//做Ifft
	void		doIfft(complex* x, int len, complex* out)
	{
		memset(&out[0], 0, sizeof(complex)* len);
		dft_1d(len, x, out, FFTW_BACKWARD/*逆向fft*/, FFTW_ESTIMATE);
		
		//normilize
		for (size_t i = 0; i < len; i++)
		{
			out[i][0] = out[i][0] / double(len);
			out[i][1] = out[i][1] / double(len);
		}
	}

	//复数相乘
	void	mult_complex(
		complex* in1, complex* in2,
		int count,
		complex* out)
	{
		for (size_t i = 0; i < count; i++)
		{
			double _Rightreal = (double)in2[i][0];
			double _Rightimag = (double)in2[i][1];

			out[i][0] = in1[i][0] * _Rightreal
				- in1[i][1] * _Rightimag;

			out[i][1] = in1[i][0] * _Rightimag
				+ in1[i][1] * _Rightreal;
		}
	}

public:

	complex*		m_fftIn;
	complex*		m_fftH;		//冲击函数fft缓存
	complex*		m_fftS;		//分段数据fft缓存
	complex*		m_fftM;		//中间相乘计算结果
	complex*		m_endC;		//缓存最后剩余的数据
	complex*		m_hC;		//缓存h
};
